{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "fdb5aa69",
   "metadata": {
    "toc": true
   },
   "source": [
    "<h1>Table of Contents<span class=\"tocSkip\"></span></h1>\n",
    "<div class=\"toc\"><ul class=\"toc-item\"><li><span><a href=\"#吉布斯采样过程\" data-toc-modified-id=\"吉布斯采样过程-1\"><span class=\"toc-item-num\">1&nbsp;&nbsp;</span>吉布斯采样过程</a></span></li><li><span><a href=\"#指标筛选过程\" data-toc-modified-id=\"指标筛选过程-2\"><span class=\"toc-item-num\">2&nbsp;&nbsp;</span>指标筛选过程</a></span></li><li><span><a href=\"#评价指标对比\" data-toc-modified-id=\"评价指标对比-3\"><span class=\"toc-item-num\">3&nbsp;&nbsp;</span>评价指标对比</a></span></li></ul></div>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e88f8b75",
   "metadata": {},
   "source": [
    "# 吉布斯采样过程"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "64001c6e",
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import shap\n",
    "import random\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import datetime as dt\n",
    "import xgboost as xgb\n",
    "import warnings\n",
    "warnings.filterwarnings(\"ignore\")\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib\n",
    "matplotlib.rcParams['font.sans-serif'] = ['SimHei']\n",
    "matplotlib.rcParams['axes.unicode_minus'] = False\n",
    "import xgboost\n",
    "from xgboost import XGBClassifier\n",
    "from sklearn.metrics import roc_curve\n",
    "from sklearn.metrics import roc_auc_score\n",
    "from sklearn.metrics import accuracy_score\n",
    "from sklearn.metrics import cohen_kappa_score\n",
    "from sklearn.model_selection import GridSearchCV\n",
    "from sklearn.model_selection import train_test_split\n",
    "from sklearn.metrics import roc_curve,accuracy_score"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ffe11682",
   "metadata": {},
   "outputs": [],
   "source": [
    "def gibbs_sampling(data,sample,turns):\n",
    "    def ic(data,sample):\n",
    "        model = xgb.XGBClassifier()\n",
    "        evalset = [(data,sample)]\n",
    "        model.fit(data,sample, eval_metric='logloss',eval_set=evalset)\n",
    "        results = model.evals_result()\n",
    "        return results['validation_0']['logloss'][-1]\n",
    "\n",
    "    # ------------------------- \n",
    "    record = {}\n",
    "    features = data.columns.tolist()\n",
    "    result = {}\n",
    "    criterion_value = {}\n",
    "    \n",
    "    # -------------------------------------------------------------------------\n",
    "    # 随机选择初始点\n",
    "    # 在实际使用中也可以从单变量模型开始\n",
    "    init = pd.Series(np.random.rand(len(features)), index = features)\n",
    "    init = init.where(init >= 15 / len(features), 1).where(init < 15 / len(features), 0).to_dict()\n",
    "    result[0] = init\n",
    "    for i in range(1, turns + 1):\n",
    "        result[i] = result[i - 1].copy()\n",
    "        record[i] = {}\n",
    "        criterion_value[i] = {}\n",
    "        for j in random.sample(features, len(features)):\n",
    "            k = 1\n",
    "            if result[i][j] == 1:\n",
    "                feature_in = pd.Series(result[i])\n",
    "\n",
    "                feature_out = result[i].copy()\n",
    "                feature_out[j] = 0\n",
    "                feature_out = pd.Series(feature_out)\n",
    "\n",
    "            else:\n",
    "                feature_out = pd.Series(result[i])\n",
    "\n",
    "                feature_in = result[i].copy()\n",
    "                feature_in[j] = 1\n",
    "                feature_in = pd.Series(feature_in)\n",
    "            \n",
    "            # -----------------------------------------------------------------\n",
    "            feature_in = feature_in[feature_in == 1].index.tolist()\n",
    "            feature_out = feature_out[feature_out == 1].index.tolist()\n",
    "            \n",
    "            if len(feature_out) > 0:\n",
    "                try:\n",
    "                # 实际使用中可能会因为进入模型的变量过多产生奇异矩阵，因此使用try来过滤这种特殊情况\n",
    "                    bic_out = ic(data.loc[:,feature_out],sample)\n",
    "                    try:\n",
    "                        bic_in = ic(data.loc[:,feature_in],sample)\n",
    "                        \n",
    "                        p = 1 / (1 + np.exp(bic_in - bic_out))\n",
    "                        \n",
    "                        rand = random.random()\n",
    "                        if rand < p:\n",
    "                            result[i][j] = 1\n",
    "                            criterion_value[i][j] = bic_in\n",
    "                            \n",
    "                        else:\n",
    "                            result[i][j] = 0\n",
    "                            criterion_value[i][j] = bic_out\n",
    "                            \n",
    "                            record[i][j] = 1\n",
    "                        print(str(100 * i / turns) + '%' + '_succeed')\n",
    "                    except:\n",
    "                        criterion_value[i][j] = bic_out\n",
    "                        result[i][j] = 0\n",
    "                        record[i][j] = 0\n",
    "                        print(str(100 * i / turns) + '%' + '_failed')\n",
    "\n",
    "                except:\n",
    "                    try:\n",
    "                        bic_in = ic(data.loc[:,feature_in],sample)\n",
    "                        criterion_value[i][j] = bic_in\n",
    "                        result[i][j] = 1\n",
    "                        record[i][j] = 1\n",
    "                    except:\n",
    "                        result[i][j] = 0\n",
    "                        record[i][j] = 0\n",
    "                        continue\n",
    "                    print(str(100 * i / turns) + '%' + '_failed')\n",
    "                    \n",
    "            else:\n",
    "                # 在剔除该变量后已经无候选变量时，\n",
    "                # 在迭代可收敛，能计算得到 bic 的情况下直接记录该特征变量进入模型\n",
    "                # 若加入该变量以后迭代不可收敛， 则该特征略过\n",
    "                try:\n",
    "                    bic_in = ic(data.loc[:,feature_in],sample)\n",
    "                    criterion_value[i][j] = bic_in\n",
    "                    result[i][j] = 1\n",
    "                    record[i][j] = 1\n",
    "                    print(str(100 * i / turns) + '%' + '_succeed')\n",
    "                except:\n",
    "                    result[i][j] = 0\n",
    "                    record[i][j] = 0\n",
    "                    print(str(100 * i / turns) + '%' + '_failed')\n",
    "                \n",
    "            k += 1\n",
    "        print(str(100 * i / turns) + '%')\n",
    "        \n",
    "    result = pd.DataFrame(result)\n",
    "    criterion_value = pd.DataFrame(criterion_value)\n",
    "    record = pd.DataFrame(record)\n",
    "    confidence = result.iloc[:,1:].sum(axis = 1) / turns\n",
    "    return [confidence, result, criterion_value, record]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0948fc94",
   "metadata": {},
   "outputs": [],
   "source": [
    "data = pd.read_csv(\"data.csv\")\n",
    "df = data.sample(n=2000, frac=None, replace=False, weights=None, random_state=0, axis=0)\n",
    "data = df.iloc[:,2:]\n",
    "data.index = list(range(data.shape[0]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a2cd485f",
   "metadata": {},
   "outputs": [],
   "source": [
    "sample = df.iloc[:,1]\n",
    "sample.index = list(range(sample.shape[0]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c284d588",
   "metadata": {},
   "outputs": [],
   "source": [
    "turns = 400 # 模拟轮数\n",
    "confidence, result, criterion_value, record = gibbs_sampling(data,sample,turns)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "126f0ca0",
   "metadata": {},
   "outputs": [],
   "source": [
    "result.to_csv(\"result_xgb.csv\",index=True)\n",
    "record.to_csv(\"record_xgb.csv\",index=True)\n",
    "criterion_value.to_csv(\"criterion_value_xgb.csv\",index=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "45022fb5",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-09-28T01:06:16.161334Z",
     "start_time": "2021-09-28T01:06:16.152960Z"
    }
   },
   "source": [
    "# 指标筛选过程"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9ce95ba5",
   "metadata": {},
   "outputs": [],
   "source": [
    "num_list = []\n",
    "index = []\n",
    "for i in np.linspace(0,1):\n",
    "    number = len(confidence[confidence>=i])\n",
    "    index.append(i)\n",
    "    num_list.append(number)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "72ca3e19",
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "plt.figure(figsize=(5,5))\n",
    "plt.plot(index,num_list,label='特征数量',color='b')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "76e4fc2e",
   "metadata": {},
   "outputs": [],
   "source": [
    "feature = list(confidence[confidence>=0.52].index)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "526349c0",
   "metadata": {},
   "outputs": [],
   "source": [
    "for i in range(len(feature)):\n",
    "    feature[i] = feature[i][:-1]\n",
    "\n",
    "feature = list(set(feature))\n",
    "len(feature)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a209a710",
   "metadata": {},
   "source": [
    "# 评价指标对比"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0ac68a8b",
   "metadata": {},
   "outputs": [],
   "source": [
    "data = pd.read_excel(\"分箱18.xlsx\",header=0)\n",
    "data.to_csv(\"分箱18.csv\",index=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "99dda45f",
   "metadata": {},
   "outputs": [],
   "source": [
    "X_before = data.iloc[:,:-1]\n",
    "Y_before = data.iloc[:,-1]\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7493bd80",
   "metadata": {},
   "outputs": [],
   "source": [
    "from collections import Counter\n",
    "print(Counter(Y_before))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "887b585e",
   "metadata": {},
   "outputs": [],
   "source": [
    "from imblearn.over_sampling import SMOTE\n",
    "smo = SMOTE(sampling_strategy={0: 261784,1:65446},random_state=0)\n",
    "X_smo, y_smo = smo.fit_resample(X_before, Y_before)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c70b3482",
   "metadata": {},
   "outputs": [],
   "source": [
    "X_smo.to_csv(\"xgb_smo_x.csv\",index=True)\n",
    "y_smo.to_csv(\"xgb_smo_y.csv\",index=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "424836d0",
   "metadata": {},
   "outputs": [],
   "source": [
    "# 设置测试集占总样本比例为0.2\n",
    "X_train_before, X_test_before, y_train_before, y_test_before = train_test_split(X_smo.astype(int), y_smo.astype(int), random_state = 0, test_size = 0.3)\n",
    "\n",
    "unique_, count_ = np.unique(y_train_before, return_counts=True)\n",
    "for unique, count in zip(unique_, count_):\n",
    "    print(\"训练集中类别为'{0}'的样本有{1}个\".format(unique, count))\n",
    "\n",
    "print()\n",
    "unique_, count_ = np.unique(y_test_before, return_counts=True)\n",
    "for unique, count in zip(unique_, count_):\n",
    "    print(\"测试集中类别为'{0}'的样本有{1}个\".format(unique, count))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4ec1fbd1",
   "metadata": {},
   "outputs": [],
   "source": [
    "xgboost_model = XGBClassifier(random_state=0, n_jobs=-1,\n",
    "                                    n_estimators=300,\n",
    "                                    max_depth=5,\n",
    "                                    subsample=1, colsample_bytree=1,\n",
    "                                    learning_rate=0.2,\n",
    "                                    scale_pos_weight=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "857b8820",
   "metadata": {},
   "outputs": [],
   "source": [
    "xgboost_model_fit = xgboost_model.fit(X_train_before, y_train_before)\n",
    "y_pred_before = xgboost_model_fit.predict(X_test_before)\n",
    "y_proba_before = xgboost_model_fit.predict_proba(X_test_before)\n",
    "acc_xgboost = accuracy_score(y_test_before,y_pred_before)\n",
    "kappa_xgboost = cohen_kappa_score(y_test_before,y_pred_before)\n",
    "auc_xgboost = roc_auc_score(y_test_before, y_proba_before[:,1])\n",
    "fpr_xgboost, tpr_xgboost, thresholds_roc_xgboost= roc_curve(y_test_before,y_proba_before[:,1])\n",
    "ks_xgboost = max([trp_ - fpr_ for trp_, fpr_ in zip(tpr_xgboost, fpr_xgboost)])\n",
    "results1 = pd.DataFrame()\n",
    "results1['准确率_1'] = [acc_xgboost]\n",
    "results1['AUC_1'] = [auc_xgboost]\n",
    "results1['KS值_1'] = [ks_xgboost]\n",
    "results1['kappa_1'] = [kappa_xgboost]\n",
    "results1.index = ['XGBoost']\n",
    "results1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e8bee27c",
   "metadata": {},
   "outputs": [],
   "source": [
    "data = pd.read_excel(\"特征30.xlsx\",header=0)\n",
    "data.to_csv(\"特征30.csv\",index=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7552d321",
   "metadata": {},
   "outputs": [],
   "source": [
    "X = data.iloc[:,:-1]\n",
    "Y = data.iloc[:,-1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9da10ac3",
   "metadata": {},
   "outputs": [],
   "source": [
    "from collections import Counter\n",
    "print(Counter(Y))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ebe42f77",
   "metadata": {},
   "outputs": [],
   "source": [
    "from imblearn.over_sampling import SMOTE\n",
    "smo = SMOTE(sampling_strategy={0: 261784,1:65446},random_state=0)\n",
    "X, y = smo.fit_resample(X, Y)\n",
    "print(Counter(y))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0b0ba824",
   "metadata": {},
   "outputs": [],
   "source": [
    "X.to_csv(\"xgb_smo_x30.csv\",index=True)\n",
    "y.to_csv(\"xgb_smo_y30.csv\",index=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "80fed8c9",
   "metadata": {},
   "outputs": [],
   "source": [
    "df = X\n",
    "X = df.loc[:,feature]\n",
    "X.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f0b51cd5",
   "metadata": {},
   "outputs": [],
   "source": [
    "X_train, X_test, y_train, y_test = train_test_split(X,y,random_state = 420,test_size = 0.3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fd369556",
   "metadata": {},
   "outputs": [],
   "source": [
    "xgboost_model = XGBClassifier(base_score=0.5,\n",
    "                        booster='gbtree',\n",
    "                        colsample_bylevel=0.8,\n",
    "                        colsample_bynode=1,\n",
    "                        colsample_bytree=0.8,\n",
    "                        gamma=0,\n",
    "                        gpu_id=-1,\n",
    "                        importance_type='gain',\n",
    "                        interaction_constraints='',\n",
    "                        learning_rate=0.1,\n",
    "                        max_delta_step=0,\n",
    "                        max_depth=3,\n",
    "                        min_child_weight=0.9,\n",
    "                        missing=np.nan,\n",
    "                        monotone_constraints='()',\n",
    "                        n_estimators=200,\n",
    "                        n_jobs=-2,\n",
    "                        num_parallel_tree=1,\n",
    "                        objective='binary:logistic',\n",
    "                        random_state=0,\n",
    "                        reg_alpha=0,\n",
    "                        reg_lambda=1,\n",
    "                        scale_pos_weight=1,\n",
    "                        subsample=0.8,\n",
    "                        tree_method='exact',\n",
    "                        validate_parameters=1,\n",
    "                        verbosity=None)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1e62a7d7",
   "metadata": {},
   "outputs": [],
   "source": [
    "xgboost_model_fit = xgboost_model.fit(X_train, y_train)\n",
    "y_pred_1 = xgboost_model_fit.predict(X_test)\n",
    "y_proba_1 = xgboost_model_fit.predict_proba(X_test)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "da6896d3",
   "metadata": {},
   "outputs": [],
   "source": [
    "acc_xgboost = accuracy_score(y_test,y_pred_1)\n",
    "kappa_xgboost = cohen_kappa_score(y_test,y_pred_1)\n",
    "auc_xgboost = roc_auc_score(y_test, y_proba_1[:,1])\n",
    "fpr_xgboost, tpr_xgboost, thresholds_roc_xgboost= roc_curve(y_test,y_proba_1[:,1])\n",
    "ks_xgboost = max([trp_ - fpr_ for trp_, fpr_ in zip(tpr_xgboost, fpr_xgboost)])\n",
    "results1 = pd.DataFrame()\n",
    "results1['准确率_1'] = [acc_xgboost]\n",
    "results1['AUC_1'] = [auc_xgboost]\n",
    "results1['KS值_1'] = [ks_xgboost]\n",
    "results1['kappa_1'] = [kappa_xgboost]\n",
    "results1.index = ['XGBoost']\n",
    "results1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0ef082db",
   "metadata": {},
   "outputs": [],
   "source": [
    "feature = [\n",
    "    'home_ownership', 'emp_length', 'dti', 'tot_cur_bal', 'num_actv_rev_tl',\n",
    "    'total_rev_hi_lim', 'verification_status', 'mo_sin_old_rev_tl_op',\n",
    "    'mort_acc', 'term', 'installment', 'num_il_tl', 'open_rv_12m'\n",
    "]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "64b25a03",
   "metadata": {},
   "outputs": [],
   "source": [
    "X = df.loc[:,feature]\n",
    "X.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e086ae45",
   "metadata": {},
   "outputs": [],
   "source": [
    "X_train, X_test, y_train, y_test = train_test_split(X,y,random_state = 420,test_size = 0.3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "40303e3a",
   "metadata": {},
   "outputs": [],
   "source": [
    "xgboost_model = XGBClassifier(base_score=0.5,\n",
    "                        booster='gbtree',\n",
    "                        colsample_bylevel=0.8,\n",
    "                        colsample_bynode=1,\n",
    "                        colsample_bytree=0.8,\n",
    "                        gamma=0,\n",
    "                        gpu_id=-1,\n",
    "                        importance_type='gain',\n",
    "                        interaction_constraints='',\n",
    "                        learning_rate=0.1,\n",
    "                        max_delta_step=0,\n",
    "                        max_depth=3,\n",
    "                        min_child_weight=0.9,\n",
    "                        missing=np.nan,\n",
    "                        monotone_constraints='()',\n",
    "                        n_estimators=200,\n",
    "                        n_jobs=-2,\n",
    "                        num_parallel_tree=1,\n",
    "                        objective='binary:logistic',\n",
    "                        random_state=0,\n",
    "                        reg_alpha=0,\n",
    "                        reg_lambda=1,\n",
    "                        scale_pos_weight=1,\n",
    "                        subsample=0.8,\n",
    "                        tree_method='exact',\n",
    "                        validate_parameters=1,\n",
    "                        verbosity=None)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1ae3c68b",
   "metadata": {},
   "outputs": [],
   "source": [
    "xgboost_model_fit = xgboost_model.fit(X_train, y_train)\n",
    "y_pred_1 = xgboost_model_fit.predict(X_test)\n",
    "y_proba_1 = xgboost_model_fit.predict_proba(X_test)\n",
    "acc_xgboost = accuracy_score(y_test,y_pred_1)\n",
    "kappa_xgboost = cohen_kappa_score(y_test,y_pred_1)\n",
    "auc_xgboost = roc_auc_score(y_test, y_proba_1[:,1])\n",
    "fpr_xgboost, tpr_xgboost, thresholds_roc_xgboost= roc_curve(y_test,y_proba_1[:,1])\n",
    "ks_xgboost = max([trp_ - fpr_ for trp_, fpr_ in zip(tpr_xgboost, fpr_xgboost)])\n",
    "results1 = pd.DataFrame()\n",
    "results1['准确率_1'] = [acc_xgboost]\n",
    "results1['AUC_1'] = [auc_xgboost]\n",
    "results1['KS值_1'] = [ks_xgboost]\n",
    "results1['kappa_1'] = [kappa_xgboost]\n",
    "results1.index = ['XGBoost']\n",
    "results1"
   ]
  }
 ],
 "metadata": {
  "hide_input": false,
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.8"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": false,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": true,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
